Scaling behavior of a square-lattice Ising model with competing interactions in a 

uniform field 
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Transfer-matrix methods, with the help of finite-size scaling and conformal invariance concepts, 
are used to investigate the critical behavior of two-dimensional square- lattice Ising spin-1/2 systems 
with first- and second-neighbor interactions, both antiferromagnetic, in a uniform external field. On 
the critical curve separating coUinearly-ordered and paramagnetic phases, our estimates of the con- 
formal anomaly c are very close to unity, indicating the presence of continuously-varying exponents. 
This is confirmed by direct calculations, which also lend support to a weak-universality picture; 
however, small but consistent deviations from the Ising-like values 77 — 1/4, 7/1^ = 7/4, P/v = 1/8 
are found. For higher fields, on the line separating row-shifted (2 x 2) and disordered phases, we 
find values of the exponent rj very close to zero. 
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I. INTRODUCTION 

The study of frustration in magnetism has been a 
very active field of research in the recent past, both 
theoretically and experimentally. While experimentally- 
realizable frustrated magnets typically have a closer cor- 
respondence to quantum (i.e., Heisenberg or XY) spin 
models than to classical, Ising-like ones, their behavior 
turns out to be rather intricate. Thus, theoretical and/or 
numerical investigation of frustrated classical spin sys- 
tems may, by virtue of their simplified character, help 
unravel some basic features which are common to frus- 
trated magnets in general. 

In this paper we investigate two-dimensional spin-1/2 
Ising systems on a square lattice with first- and second- 
neighbor couplings, both antiferromagnetic, in the pres- 
ence of a uniform magnetic field. The Hamiltonian is 
given by: 



n 



•h 2_^ cTi aj + J2 2_^ ^i ^3 



HY^a, 



(1) 



NN 
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where Ji, J2 > 0, NN and NNN stand respectively for 
next-neighbor and next-nearest-neighbor pairs, and ai — 
±1. Here, all fields, coupling strengths and temperatures 
are given in units of Ji, unless otherwise stated. We have 
kept J2 = 1 in all calculations reported in this work, 
except at the end of Sec. IIII Al where J2 = 2 and 0.75 
were briefly considered (both for H = 0). 

In line with the initial considerations given above, this 
may be considered a classical approximation for the Ji — 
J2 (Heisenberg) model [1, 2]. Nevertheless, as shown in 
the following, the model described by Eq. ([T]) exhibits 
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intricate features of its own, several of which are not fully 
understood so far. Depending oir the relative streirgth 
of the associated parameters, such setup of competing 
interactions can generate various types of ordered phases 
at low temperature; more often than not, the transitions 
between these and the high-temperature (paramagnetic) 
state do not belong to the standard Ising universality 
class 0,[1|. 

The problem studied here has been analyzed by several 
numerical techniques in the past; Refs. i3, and J, provide 
excellent summaries of earlier work, as well as illustra- 
tions of the use of up-to-date Monte-Carlo (MC) simula- 
tion techniques for this case. 

We use transfer- matrix (TM) methods [5|, in con- 
junction with finite-size scaling ^] and conformal in- 
variance [3| concepts, to determine the location of the 
phase boundaries of systems described by Eq. ([T]), and 
the universality classes of the associated phase transi- 
tions. TM methods, especially in the strip geometry used 
in this work, are to some extent complementary to MC 
simulations, in that they provide straightforward proce- 
dures for evaluation of the conformal anomaly, or central 
charge [8|, as well as the decay-of-correlations exponent 
r] (via the amplitude-exponent relationship [9]). Both 
quantities play an important role in the ideirtification of 
the universality classes pertaining to phase trairsitions 
in two-dimensional systems, and neither is directly ac- 
cessible via MC methods (although direct estimates of 
rj can be produced by following the decay of spin-spin 
correlations with distance in an MC context, no simple 
relationship applies, such as the one given by conformal 
invariance on strips |9|]). On the other hand, similarly to 
MC techniques, TM calculations also provide estimates 
of critical temperatures, specific heats, magnetizations, 
and susceptibilities. 

Iir Section [H] we recall the calculational methods used, 
as well as the finite-size scaling concepts and techniques 



employed in the analysis of our results. Our numerical 
results ioT H = are given in Section UlIAi and those for 
H y^ in Section IIIIBI Finally, in Sec. llVi concluding 
remarks are made. 



II. CALCULATIONAL METHOD AND 
FINITE-SIZE SCALING 

We set up the TM on strips of width N sites, with pe- 
riodic boundary conditions across. The coordinate axes 
coincide with the directions of the first-neighbor bonds. 
We used 4 < A*" < 22. For comparison, earlier TM stud- 
ies of this problem [l^, [ll| could only reach iV = 12. 
For the case J2 = 1/4 in zero external field, results for 
iV < 18 are available ^. 

With Ai, A2 being the two largest eigenvalues (in ab- 
solute value) of the TM, the dimensionless free energy 
per site is given by /Ar(T) = iV~^lnAi, while kn{T) — 
In I A1/A2 I is the inverse correlation length on a strip of 
width N sites |5l]. 

It must be stressed that we do not make any assump- 
tions about symmetry properties of the TM's eigenvec- 
tors. Starting from the full set of 2^ basis vectors, 
the eigenvector | ui), corresponding to Ai, is isolated by 
the power method, while | V2) is again extracted via the 
power method, combined with repeated Gram-Schmidt 
orthogonalization to \vi). This way we make sure that 
the most strongly diverging correlation length is evalu- 
ated, that is, the one which truly corresponds to the order 
parameter for the transition under scrutiny. This is es- 
pecially relevant in the present case, where critical lines 
corresponding to order parameters of differing symme- 
tries can become very close (see Section PlI Bl below) . 

Here we assume that the transition is always of second 
order, which is implicit in the statements made in the 
preceding paragraph. By now this seems out of doubt, 
at least for j/2 > 1 which is the parameter range of in- 
terest here [3, |j|- For fixed J2 and H, say, we locate 
the approximate (A'^-dependent) critical temperature T^ 
by solving the basic equation of the phenomenological 
renormalization group (PRG) Q: 



Ni^n{T)=N'kn'{T) 



(2) 



Depending on the shape of the critical curve, it may 
be more convenient to keep T fixed and vary H, in 
which case the critical behavior is expressed in terms of 
\H-Hc\,a.nd Eq. ^ gives H^j. The strip widths N and 
N' are to be taken as close as possible for improved con- 
vergence of results against increasing N. In order to obey 
ground-state symmetries, here we used only even N, N', 
so N' = N - 2. For H = and J2 < 1/2, in which case 
the ordered phase is Neel-like ^], one can use both odd 
and even TV; indeed, good results are found from PRG 
with N' = N - 1 in that region [T2|. For J2 > 1/2, we 
found that: (i) the latter procedure does not give physi- 
cally meaningful solutions for Eq. ([2]); and (ii) although 
PRG with both N and N' odd gives the same limiting Tc 



with iV ^- cxD as when both strip widths are even (albeit 
with much slower convergence), estimates of quantities 
other than the critical temperature are unreliable. 

Estimates of the thermal exponent i/t = ^/v are given 
by J5i]: 



yr 



1 



\n{N/N') 



(3) 



where k^, k^, are temperature derivatives of the inverse 
correlation lengths, taken at T^. Finite- A^ estimates of 
the exponent rj are given by the conformal invariance 
relation ^9||: 



1]N = TT ^7VKAr(r^) 



(4) 



The convergence of finite-TV approximants given by 
Eqs. (12), (131), and ^ towards their iV — > 00 values has 
been extensively discussed (l3l - [l7l |. For (unfrustrated) 
Ising-like systems on strips with periodic boundary con- 
ditions across, the rate of convergence goes like 



X 



N 



x^ 



aN- 



(5) 



with w « 3 for X = T* , and a; « 2 for X — yx (for some 
simple cases, this can be shown analytically [13, 16, 17]). 
By taking sets of three successive finite- iV estimates, one 
can use w as an adjustable parameter in Eq. ([S]), and 
produce a new, shorter, sequence which can then be it- 
erated again, and so on. Such iterated three-point fit 
technique can produce very accurate final estimates of 
critical quantities [Tj, [l^, \19^ . 

Once Tc is found, as described above, to good accu- 
racy (or if its exact value is known, for example via du- 
ality arguments ^^), sequences of assorted quantities 
can be evaluated at the extrapolated critical point, for 
increasing iV. From these, one can usually extract es- 
timates of critical exponents which converge faster and 
more smoothly than if the calculations were done at the 
respective pseudo-critical temperatures T^ |13i, ,2(| . One 
is interested in (per site) specific heats, susceptibilities, 
and magnetizations, which behave as [6]: 






(6) 



Both Cjv and xn are found from suitable second deriva- 
tives of the free energy [18|. The exponent ratio a/i^ can 
then be extracted from three-point fits of sequences of 
Cn, as explained in connection with Eq. ([5]). For ^/ly, 
one initially obtains a sequence of exponent estimates via 
two-point fits of susceptibility data, and then proceeds to 
extrapolating such a sequence via three-point fits \X^- 

The spontaneous magnetization rriN is difficult to cal- 
culate in a finite-size scaling context, because it is iden- 
tically zero for a finite system. For quan tum chains at 
T = this problem can be overcome |2l| . by exploiting 
the fact that there the largest eigenvalue of the TM gives 



the internal energy: in a first-order degenerate perturba- 
tion scheme, appropriate consideration of non-diagonal 
matrix elements enables one to extract the magnetiza- 
tion in the zero-field limit. For classical spins on strips, 
the corresponding eigenvalue of the TM gives the free 
energy instead, and the perturbation-theory procedure 
used for quantum systems [2ll] cannot be translated to 
our case. 

In this work we estimated the finite-size magnetiza- 
tion exponent p/v by calculating the average squared 
magnetization per column at Tc, (M^). Considering, for 
example, a ferromagnet, denoting by cr = {cri • • • aj^} 
the 2^^ column basis states, and with '0(c)j i^{^) being 
respectively the dominant left and right eigenvectors of 
the TM, one has [i|: 



(M2) = 



T.J{^){T.l,^r)\{'^) 



EX^)V'(a) 



At the critical point, one should have: 



1 

TV 



(A/2)l/2_^-,9/. 



(7) 



(8) 



For a square-lattice antiferromagnet with only first- 
neighbor interactions, Ui in Eq. ^ must be replaced by 
(— l)'cri, so the staggered character of the order param- 
eter is properly taken into account. The corresponding 
adaptation for the system of interest here is discussed 
in Section IIIII below. We tested this procedure, with 
the appropriate (uniform or staggered) version of the 
magnetization, on the following square-lattice Ising sys- 
tems: (i) ferromagnet with nearest-neighbor couplings 
only; (ii) ferromagnet with first- and second-neighbor in- 
teractions, J2 — 1] and (iii) antiferromagnet with first- 
and second-neighbor interactions, J2 — 1/4 (for which 
the ordered state is Neel-like [l^l)- In all three cases, T^ 
is known either exactly or to a very good approximation, 
and the transition is in the Ising universality class [12], so 
j3/v — 1/8. In order to gauge the likely systematic errors 
for our intended final application (see Section |III[) . we 
considered 4 < iV < 22, and only even N . All resulting 
sequences gave estimates of P/v monotonically growing 
with iV, pointing to extrapolated values between 0.1240 
and 0.1247, so the systematic error is less than 1% for 
this range of TV. 

Another quantity of interest to be calculated at Tc is 
the conformal anomaly c, given by the N^'^ finite-size 
correction of the critical free energy per site [8]. For the 
present case of strips with periodic boundary conditions 
across, one has: 



the model studied here belongs to the latter category (for 
H = 0, this has already been pointed out in Ref. |24[ ). 

Additionally, one can both double-check the robustness 
of extrapolations of Tc and yr from Eqs. ^ and ([3]), as 
well as investigate the other quantities of interest, by 
scanning the neighborhood of the critical point with the 
help of finite-size scaling ideas [6]. Taking, for instance, 
77jv(r) = 7r~^ NknJT), and allowing for corrections to 
scaling, we write |25l [2g : 

m{T) = f{u) + N-'^g{u) , u = 7Vi/'^(T - Tc) . (10) 

where a; > is the exponent associated with the leading 
irrelevant operator [ see Eq. ^ ] . Close enough to Tc 
the scaling functions in Eq. ([TU)) should be amenable to 
Taylor expansions. One has: 






N- 



A:=0 



6fcW 



(11) 



where 77 is to be compared with the N —> 00 extrapolated 
value of the 77 jv of Eq. ([4]) . 

One looks for values of T^, ^, w and the {aj, bk} which 
optimize data collapse upon plotting rj^iT) — N~'^g{u) 
against u. In practice, good fits are generally found with 
jm, km not exceeding 2 or 3 |25l.l26l|. 

Considering now the finite-size susceptibility xn{T)^ 
finite-size scaling Q suggests a form 

XN{T) = N^'''fA^), (12) 

where 7 is the susceptibility exponent. Following Refs.l25l 
and 1261 we write (again, allowing for corrections to scal- 
ing): 



lnx7v(T) = -lniV + Va,^u^+iV- 

^-1 



k=0 



(13) 



In order to reduce the number of fitting parameters, it is 
usual to keep 1/v and Tc fixed at their central estimates 
obtained, e.g., via Eqs. (fTU)) and (fTTj) . allowing ^/v to 
vary. 

Expressions similar to Eq. (J13p can be written for mag- 
netizations and specific heats, yielding estimates of the 
exponents 13 /v and a/ v. 

In the present context, one should interpret the expo- 
nent u) in Eqs. (O, pi?)) , and (|13p as an effective one, 
representing all orders of corrections to scaling (which 
may also turn out to have rather different amplitudes 
for different quantities). Thus, in practice a somewhat 
broad range of results (say, 1 < 1^ ;$ 3) can be accepted 
when considering data collapse optimization for distinct 
quantities related to the same problem. 



lM{Tc) = h + ^,+0{N-') 



(9) 



While models with c < 1 are associated with universality 
classes with fixed values of the critical exponents, those 
with c > 1 can have continuously varying exponents [22|, 
[23fl . As shown below, there are strong indications that 



III. NUMERICAL RESULTS 

A. H = 

In zero external field, for J2 < 1/2 the ground-state or- 
dering is of the Neel type, with the two sublattices aligned 



antiparallel to each other; the transition is second-order, 
in the Ising universahty class [12|- At J2 = 1/2 the 
ground state is macroscopically degenerate, and the crit- 
ical temperature is zero \4\. For J2 > 1/2 the lowest 
energy corresponds to collinear order, with alternating 
rows (or columns) of parallel spins. For J2 > 1/2 the 
transition is first order, and evidence has been found 
that it remains so, at least up to J2 ~ 0.9 [4^, i2J|. As 
J2 increases further, the second-order character returns. 
The bulk of extant evidence [1, [H ^^ indicates that 
the transition is second order for J2 = 1. Finally, for 
J2 S> 1 one has a picture of two weakly-coupled antifer- 
romagnetic lattices, thus in this limit Tc/ J2 approaches 
the Ising value, 2/ ln(l -f %/2). 

For J2 = 1 , recent estimates of the critical temperature 
are: T^ = 2.0823(17) [l^; 2.0838(5) [lil; 2.0820(4) [3; 
and 2.0839(12) ^. By solving Eq. Q, we found a 
well-behaved sequence of T^ values, extrapolating to 
Tc = 2.08195(5) via three-point fits, with w « 3. The 
sequences for j/t and rj from Eqs. ^ and ^ extrapolate 
respectively to j/t = 1.188(2) and rj = 0.2341(1), again 
via three-point fits. 

Evaluating 7yjv(T) in the region around Tc , and em- 
ploying Eqs. (dni) and mi), gave Tc = 2.08197(5); yr = 
1.182(3), and 77 = 0.2342(1), with uj w 1.9. We used 
jm = 2, /cm = 1 in Eq. (fTTj) . Thus there is a satisfac- 
tory degree of consistency between the two methods of 
evaluation of critical quantities. 

Our numerical value for v = y^^ = 0.844(4) [from 
averaging over the two results abo ve ] is to be compared 
to iy = 0.8292(24) [27]; 0.8481(2)[2|; 0.847(4) [1|; and 
0.847(1) 129]. The value 77 = 0.20(1) was found by direct 
MC evaluation of critical correlation functions on N x N 
geometries |24| . 

By evaluating quantities at the extrapolated Tc, we 
found 77 — 0.23415(5); as expected, this is even more ac- 
curate than extrapolating the sequence of finite- A^ values 
estimated at the respective fixed points T^. 

For calculation of the zero-field susceptibility, the spe- 
cific properties of the collinear order parameter were 
taken into account as follows. For a fixed coordinate di- 
rection, say X, along which the TM proceeds, the critical 
wavevector is degenerate, being either (tt/o) x or (vr/a) y, 
with a being the lattice parameter. One thus has to take 
both (equally probable and mutually exclusive) possibil- 
ities into account and average the partial contributions 
given by each. As might be expected, we found both con- 
tributions to be of similar amplitudes (within « 10% of 
each other for fixed N); separate fits of each to power- 
law forms gave apparent exponents differing by less than 
1%. The latter discrepancy can be ascribed to residual 
lattice effects, and is expected to vanish for larger N, out 
of reach of our TM implementations at present. 

Our final result was 7/;^ — 1.772(1). Estimating ^/v 
via Eq. p^ gave 7/1^ = 1.775(1), slightly higher than the 
previous estimate but within three (rather narrow) error 
bars. Since the uncertainties quoted refer exclusively to 
the fitting procedures, i.e., no account is taken of likely 



systematic errors, one might err on the side of caution 
and allow for somewhat large uncertainties. Averaging 
over the two values found, we quote j/i' — 1.773(4). This 
way our results might be considered marginally compat- 
ible with 7/1/ — 1.750(12), quoted in Ref. 3, although it 
seems much harder to stretch our error bars to include 
the value 7/4, which would be consistent with a weak- 
universality picture of "//v, P/v, (2 — a)/^ sticking to 
the respective Ising values [30] ■ Ref. [23 quotes the range 
1.71 — 1.79 for 7/1^, based on three diferent fitting meth- 
ods. 

For the calculation of /3/i', we used as column magne- 
tization in Eq. ([7]) the following quantity: 



(M2) = i[(Af2) + (A4^] , 



(14) 



where M^, Mgt are respectively uniform and staggered 
column magnetization. This choice reflects the collinear 
nature of the ground state, with its orientational degen- 
eracy, and closely corresponds to the order parameter 
used in the MC simulations of Ref. y. See the arguments 
invoked above for the susceptibility calculation. Simi- 
larly to the test cases described in Section |lll we found 
exponent estimates monotonically growing with N; the 
extrapolated result is /3/z/ = 0.121(2), where an ad hoc 
doubling of the uncertainty found from fits has been in- 
corporated, in order to allow for the small bias shown in 
tests. 

We also evaluated critical specific heats. Finite-size 
specific-heat sequences can prove unwieldy to extrapo- 
late, even when the exponent a is positive, as in the 
case of the two-dimensional three-state Potts ferromag- 
net [ll]. Here, three-point fits of iV, iV - 2, iV - 4 data 
gave a/v increasing from « 0.31 {N = 10) to w 0.33 
{N = 20), albeit with small oscillations; a quadratic fit 
of such values against 1/A^ then gave a/i^ = 0.351(12). 

The above results are to be compared to P/v = 
0.122(4), a/iy = 0.357(8), both from Ref. "3, and a/iy = 
0.412(5) [27]. Recalling the Rushbrooke scaling relation, 
a -f 2/3 + 7 = 2, our estimates give a/v + 2(3/^ + ^jv — 
2.366(13), while 21 v = 2.370(10). Similarly, one has 
(2 — a)/v — 2.020(18). Given the rather large uncertain- 
ties found in the analysis of specific heat behavior, we do 
not believe that any actual breakdown of hyperscaling is 
present. 

We calculated free energies at T^ and fitted them to a 
quadratic form in l/N"^ , thus extracting estimates of the 
conformal anomaly [a]. From fits of data in the range 
[A',„in,22], with 4 < A^min < 14, we found estimates of 
c decreasing monotonically from 1.074(1) for A'min = 4 
to 1.056(1) for A^min = 14. Uncertainties quoted relate 
exclusively to the fitting procedure. This range of esti- 
mates cornpares favorably with the corresponding result 
from Ref. [2J, c = 1.0613(6). It must be kept in mind 
that what one is seeing most likely amounts to strong 
crossover effects distorting a picture where c = 1 |24| . 

We also considered J2 = 2, for which case we obtained, 
from extrapolating sequences generated via Eqs. ©, 



dSD, and (H , TJJ2 = 2.2248(1), yr = 1.052(1), rj = 
0.2391(3). Recent results are: Tc/J2 = 2.226(5) [23], and 
2.227(5) (Hi; yr = 1.066(1) IHl. Evaluation of finite-size 
susceptibilities, magnetizations, and specific heats at Tc 
gave 7/1/ = 1.756(2); P/iy = 0.120(2); a/v < 0.1 (esti- 
mates for the latter quantity were plagued by the same 
sort of irregularities reported for J2 — 1 above). Overall, 
these values are consistent with a picture of continuously- 
varying ex pon ents, approaching the Ising ones as J2 in- 
creases [^, l2J|. Conformal-anomaly estimates are very 
close to c = 1.010; again, this is consistent with the trend 
towards c = 1, followed by fitted results upon increasing 
J2, found in Ref. [13. 

Finally, we made J2 — 0.75, which is expected to cor- 
respond to a first-order transition [24], thus in principle 
the ideas behind Eq. ([2|) do not apply. Indeed, instead of 
varying monotonically with increasing N^ the solutions 
of Eq. © initially went up, to Tc w 1.432 at A^ = 12, 
and then became approximately constant for larger N; 
the 77 estimate from Eq. Q also initially increased, up to 
« 0.243 at TV = 14 and 16, then started decreasing for 
larger N. This indicates a correlation length which at the 
very least grows slower than N, and possibly saturates at 
scales which are out of reach of our TM calculations, that 
is, a weakly first-order transition [J, [2J]. Notwithstand- 
ing the lack of conceptual justification for using Eq. (0), 
it should be noted that Tc « 1.43 is in rather good agree- 
ment with MC estimates (see Figure 3 of Ref. |j) . 



B. H ^0 

For J2 = 1, the ground state is still collinear for iJ < 4, 
whereas for 4 < 7f < 8 it becomes a row-shifted (2 x 2) 
state 13\ . The latter consists of alternating ferro- and an- 
tiferromagnetically ordered rows (or columns), with the 
ferromagnetic ones parallel to the field. The added de- 
gree of freedom (relative to a 2 x 2 state) is that the 
antiferromagnetic chains can slide freely relative to each 
other, at zero energy cost. At iJ = 4 and 8, Tc = 
because of macroscopic ground-state degeneracy. The 
maximum critical temperature for the transition between 
row-shifted (2x2) order and the paramagnetic phase has 
been estimated as « 0.73, at i? « 6 [3]. 

In order to make contact with previous results, 
we initially considered two points on the collinear- 
paramagnetic transition line, respectively at H — 2.5 and 
3.3 Q. 

For H ~ 2.5, from extrapolating sequences generated 
via Eqs. dH), ©, and ^, we found Tc = 1.6846(1), 
77 — 0.2335(1). A noticeable trend reversal was ob- 
served for yr; after decreasing from 1.08 to 1.072 be- 
tween N = 8 and 12, it starts increasing smoothly, reach- 
ing 1.075 at A^ = 22. Extrapolating iV > 12 data, we 
found j/T = 1.090(4). Evaluating susceptibilities at the 
extrapolated Tc resulted in 7/1/ = 1.779(4), and anal- 
ysis of magnetizations gave (3/iy — 0.122(1). Ref. y 
gives: Tc = 1.6852(3); yr = 1.056(8); j/i^ = 1.750(14), 



P/v = 0.118(3). Finally, the conformal anomaly was es- 
timated as c = 1.066(1). 

Following the same procedure as above, we obtained 
for H = 3.3: Tc = 1.3331(5); yr = 0.940(3) (this 
time with no trend reversal upon increasing N); rj = 
0.2338(3). Finite-size susceptibility scaling at Tc gave 
7/:^ — 1.781(5), and magnetizations, (3/i' = 0.135(2). 
The evolution of /3/i> along the collinear-paramagnetic 
phase boundary is analyzed towards the end of this Sec- 
tion. Ref. i quotes: T^ = 1.3335(6); yr = 0.907(7); 
7/i^ = 1.751(14), l3/u = 0.130(5). Our estimate for the 
conformal anomaly is c = 1.042(1). 

The above results show very good numerical agreement 
with existing ones as regards critical temperatures; also, 
a picture of continuously-varying exponents such as v 
and 7 is confirmed, both directly and from the conformal 
anomaly results, which are very close to unity. Our values 
for 7/z/ and 77 do seem consistent with a weak-universality 
scenario along the collinear-paramagnetic phase bound- 
ary; however, they indicate small but consistent devia- 
tions from the Ising-like picture of 7/1^ = 7/4, r/ = 1/4. 

In order to investigate the latter point in detail, we 
proceeded to evaluating 77 and c along the full extent of 
the phase boundary. 

We first extrapolated, to iV — >■ cx), the (T, H) values ob- 
tained for sequences of solutions of Eq. ^ , with increas- 
ing N and either H fixed, or (typically for lower temper- 
atures) T fixed. This was done both for the collinear- 
paramagnetic critical line and for that separating the 
row-shifted (2 x 2) and paramagnetic phases. In the 
former case, we fitted finite- A^ data with N — 14 — 20 
to the single-power form, Eq. ([5]), finding good conver- 
gence with 3 < a; < 4 everywhere on the phase bound- 
ary. For the latter, we had extremely slow convergence 
of our PRG calculations, which limited us in practice to 
N < 16. Probably (at least partially) as a consequence 
of this, the single-power form produced rather low ad- 
justed values of w, in the neighborhood of 0.5, which 
is usually interpreted as indicating strong corrections to 
scaling. We thus resorted to ad hoc parabolic fits in N~'^, 
adjusting iV = 10 — 16 data to this latter form. 

The behavior in the low-temperature region near H = 
4, where the two distinct phase boundaries become 
close to each other, is of special interest since it has 
been sug gested that an XY-\ike region might be present 
there p3l.l3l|. 

In Figure[T]we show our results for the low-temperature 
part of the phase diagram, near H — A. Although nu- 
merical convergence difficulties prevented us from reach- 
ing T < 0.2 for the largest strip widths, we managed 
to evince clearly-defined trends followed along both crit- 
ical lines, on their approach to T = 0. Below T = 0.4, 
both the collinear-paramagnetic and row-shifted (2 x 2)- 
paramagnetic boundaries are, to a very good approxi- 
mation, straight, and pointing towards {T,H) — (0,4) 
with respective slopes 0.267(1) and 2.713(5). Concurring 
with Ref. |3|, it appears very unlikely that an XF-like 
region, or a bicritical point at T > 0, is present. In- 
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Figure 1. (Color online) Phase diagram near H = 4, show- 
ing phase boundaries: collinear-paramagnetic, with reentrant 
behavior, and row-shifted [R-s] (2 x 2)-paramagnetic. The 
points are the results of extrapolating sequences of solutions of 
Eq. (12J, obtained with fixed T and variable H (see text). Un- 
certainties are smaller than symbol sizes. Each of the dashed 
lines at r < 0.2 is the continuation of the best-fitting straight 
line joining points at 0.2 < T < 0.4, on the respective phase 
boundary. 



Figure 2. (Color online) Phase diagram for high fields 4 < 
H < 8, showing row-shifted [R-s] (2 x 2)-paramagnetic phase 
boundary: finite-TV solutions of Eq. ([2]) for A'' = 10 and 16, 
and results of extrapolation of A'^ = 10 — 16 curves (see text). 
Uncertainties in the latter are smaller than symbol sizes. Each 
of the dashed lines at T < 0.2 is the continuation of the best- 
fitting straight line joining points at 0.2 < T < 0.4, on the 
respective section of the extrapolated phase boundary. 



stead, all indications from our results are consistent with 
both critical Hues joining at a single bicritical point at 
T = 0. Additionally, we found the maximum of the reen- 
trance on the collinear-paramagnetic phase boundary to 
he H = 4.121(2) at T = 0.5. At T = 0.7 we estimate 
H = 4.060(6). These are in rather good agreement with 
the respective values H — 4.07(2) and 4.052(7), quoted 
in Ref.tl. 

Near H = 0, the phase boundary has the expected 
parabolic shape [33.l33|: 



T,{H)^T,{0)-aH^ {H ^ 0) 



(15) 



where we found a ~ 0.0595(3) by fitting data correspond- 
ing to < if < 0.4. 

As noted previously in Ref. [Tl| and above, we generally 
found finite-size effects to be much larger for the high- 
field (4 < H < 8) part of the phase diagram. This is 
illustrated in Figure[21 where finite- iV curves with the so- 
lutions of Eq. 1^ ioT N — 10 and 16 are displayed jointly 
with our final extrapolation (as described above). At 
variance with Ref. 5, where Tc = 0.7293(7) is reported at 
H = 6, our extrapolated value is Tc{H = 6) = 0.589(4). 
Although alternative procedures to our ad hoc parabolic 
extrapolations against 7V~^ can certainly be devised, it 
must be noted that for this same field intensity the so- 
lution of Eq. dH) is T^ = 0.706 aheady for N = 10, and 



decreases systematically with increasing N. 

At {T,H) = (0,8), the initial slope S = {dH^/dT)T=o 
of the critical curve gives an estimate of the reduced crit- 
ical chemical potential fi/ksTc for the hard-square lat- 
tice gas with first- and second-neighbor exclusion, via 
fi/kBTc — ~2S [11]. Our result for the chemical potential 
is 5.42, to be compared with 4.70 flT], and 4.91 [si]. This 
may indicate that our extrapolation procedures slightly 
underestimate the extent of the row-shifted (2 x 2) phase. 

Figure [3] shows our results for 77 and c along the ex- 
trapolated location of the collinear-paramagnetic border, 
parametrized by T. These are obtained from Eqs. ^ 
and ^, respectively. In part (a), comparison between 
A^ = 10 estimates and the final iV — >■ 00 extrapolation 
illustrates that residual finite-size effects contribute to- 
wards overestimating the exponent 77, for all < T < 1.78 
(the approximate point where all finite- A^ curves cross). 
On the other hand, for higher T the finite-size corrections 
change sign. The extrapolated rjxT curve is to a large ex- 
tent horizontal, near both the low- and high-temperature 
ends of the phase boundary. We estimate 77 = 0.2476(3) 
for T < 0.5, and 77 = 0.2342(3) for T > 1.3. In the inter- 
mediate region there is a crossover which becomes rather 
sharp around T = 0.8, at the upper end of the reentrant 
part of the phase diagram. The conformal-anomaly es- 
timates in part (b) show the same behavior found for 
H = in Section IIII Al and in Ref. [H, in that they are 
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Figure 3. (Color online) (a) Decay-of-correlations exponent 
ri and (b) confornial anomaly c, calculated along the extrap- 
olated coUinear-paramagnetic phase boundary. In (a), the 
dashed line gives estimates of rj from strips of width A'^ = 10 
sites, while points are extrapolations from sets of finite- A'' es- 
timates, A'^ = 14 — 20. In (b), values of c are estimated from 
quadratic fits of free-energy data against A*'"^ for TV = 14 — 20. 
Uncertainties are smaller than symbol sizes. 



Figure 4. (Color online) Finite-size magnetization exponent 
jS/f along coUinear-paramagnetic (PM) phase boundary. The 
range of fields, < H < 4.027, on the horizontal axis cor- 
responds, respectively, to 2.0820 > T > 0.75 (see text). 
Each point is the result of fitting finite-size data in the range 
12 < A/' < 20 to a single-power law, {M^)^/^ ~ N^''^/''. See 
Eqs. ® and (O. 



always slightly above unity. Similarly to the H = case, 
we also found that, upon fitting free-energy data in the 
range [A^min, 20], the estimates of c always decrease upon 
increasing iVmin- Thus, an interpretation of the present 
results as consistent with c = 1, albeit affected by strong 
crossover effects, seems credible. 

We evaluated critical magnetizations, as given by 
Eq. p^ , along the coUinear-paramagnetic phase bound- 
ary. Our calculations did not converge for T < 0.75, 
which approximately coincides with the reentrant region. 
Thus, for the part of the critical boundary where we man- 
aged to produce estimates of /3/i^, there is a one-to-one 
correspondence between field and temperature. Our re- 
sults are shown in Figure 21 parametrized by H. This 
way, it is easier to follow the evolution of quantities for 
low fields than if we used T for the horizontal axis, be- 
cause of the parabolic shape assumed by the critical curve 
in that region. One sees that the quality of fits gener- 
ally deteriorates as H increases; the shallow dip around 
H ss 1.5 is possibly related to slight inaccuracies in the 
determination of the extrapolated critical line in that re- 
gion. A more persistent trend is that towards increasing 
values for larger H. We interpret this as signalling the on- 
set of the physical effects which give rise to reentrant be- 
havior for even larger fields. Indeed, a plausible explana- 
tion for the reentrance is, to quote Ref. 3;', 'the appearance 
of (2 X 2) "clusters" that help to sustain the [collinear] 



order at low temperatures even when the external field is 
slightly bigger than 4'. This sort of cluster is not taken 
into account in our column-magnetization calculations, 
see Eq. p^ . We also know that the general effect of ne- 
glecting relevant contributions to the magnetization is to 
increase the apparent value of P/v; for instance, if M^ is 
discarded in Eq. ((M]), the estimate of /?/t/ at T = 2.0820, 
H = goes from 0.120(2) to 0.135(2). According to this 
interpretation, for H ~ 2.5 or thereabouts, (2 x 2) con- 
figurations which are locally energetically favorable start 
contributing to ordering in the iV — >■ cxd limit, but are 
not captured in the scheme of Eq. ([13]). With decreas- 
ing T and increasing H, the effect of such configurations 
becomes more relevant, providing a mechanism through 
which the apparent exponent increases, although the real 
one, we conjecture, possibly increases a little but stays 
slightly below 1/8. This would be in line with the behav- 
ior of ry, depicted in part (a) of Figure [3] 

Calculation of the thermal exponent yx via Eq. ([3]) in 
the reentrant region gave negative values, an artifact al- 
ready noticed in Refs. [l^ and|Tl|. However, evaluation of 
critical finite-size susceptibilities at T = 0.35 gives j/i' 
in the range 1.71 — 1.80, depending on the details of cor- 
rections to scaling assumed for data fitting. Although 
lacking in accuracy, this range of values is broadly con- 
sistent both with the corresponding rj estimate, and with 
the hypothesis that critical behavior obeys weak univer- 
sality all along the coUinear-paramagnetic critical line. 
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Figure 5. (Color online) Decay-of-correlations exponent 77, 
evaluated via Eq. Q , along approximate row-shifted (2 x 2)- 
paramagnetic transition lines, obtained by solving Eq. ((2]) for 
A'' as indicated, for high fields 4 < _ff < 8. 



Turning now to the high-field part of the phase dia- 
gram, we illustrate in Figure[5]our results for 77 along the 
approximate critical lines, the latter obtained by solving 
Eq. dH) for iV = 10 - 16. The curve for A^ = 12 corre- 
sponds to that shown in Figure 2 of Ref. [ill By com- 
paring the evolution of t] along the approximate critical 
lines with the evolution of the lines themselves against 
increasing N^ shown in Figure [21 one anticipates that 
calculating 77 on the extrapolated phase boundary will 
give results very close to zero, or even slightly negative. 
Indeed, this was what we found. It would appear that 
this is at least partly because our extrapolation proce- 
dures underestimate the extent of the row-shifted (2 x 2) 
phase. Evaluation of c along the extrapolated critical line 
also gave physically inconsistent results. 

Even though we are not able to produce numerically 
accurate estimates of ry for the high- field part of the phase 
diagram, the gist of the results shown in Figure [S] is that 
this must be below 0.01, and possibly even zero. We 
return to this point in the next Section. 



IV. DISCUSSION AND CONCLUSIONS 

For the model described by Eq. ([T]), with J2 = 1, 
we have established a physical picture for the coUinear- 
paramagnetic phase boundary, which is consistent with 
continuously-varying exponents along the critical line. 
Together with various pieces of numerical evidence, col- 
lected at selected points, overall support for this is given 
by the conformal anomaly results depicted in part (b) of 



Figure [3] 

There is also clear evidence that such continuously- 
varying exponents satisfy, at least approximately, a weak- 
universality scenario '30|. However, as shown for the ex- 
ponent r] in part (a) of Figure [31 our results indicate 
small but consistent deviations from the corresponding 
Ising values. 

Furthermore, such deviations are internally consistent, 
in the sense that both 77 and P/v take on values lower 
than the Ising ones, while 7/1/ is always found to be 
higher than the Ising result. For (3/h', the apparent re- 
versal of this trend found for H > 2.5 has been explained 
in Section IIIIBl as a likely effect of the same sort of lo- 
cally stable (2 x 2) configurations which, at lower T and 
higher H , become significant enough to induce reentrant 
behavior. 

Notwithstanding the compensation just referred to, 
our estimates of (7/7^) + rj in general exceed 2, though 
never by more than 2 — 3 times the respective (com- 
bined) error bar. However, it must be recalled that here 
j/iy is essentially one order of magnitude larger than rj, 
and both quantities have similar relative uncertainties, 
thus the calculated combined uncertainty is practically 
only that associated with "f/v. In such circumstances, 
the apparent violation of a fundamental scaling relation 
reflects the fact that the relative uncertainty in j/i/ was 
estimated as « 2 parts in 10^. Had this been doubled, 
all the basic conclusions from this work would still stand, 
and the mismatch would be essentially lost within the re- 
vised error bars. 

Returning to 77 as displayed in part (a) of Figure O 
the small but consistent shift between the high- and low- 
T approximately constant values [respectively, 0.2342(3) 
and 0.2476(3) ] indicates a crossover between two distinct 
weak-universality classes. Such small variations could 
probably be accounted for in the context of compacti- 
fied boson theory |23l l35|. in which continuously- varying 
critical indices are put in direct correspondence with the 
(also continuously-varying) radius R associated with the 
underlying field theory. 

In general, both for iJ = and H ^ Q (the latter, 
along the collinear-paramagnetic critical line) our results 
for the location of critical points, and exponents such 
as yT^ l/v, and 13/v, are mostly compatible, within er- 
ror bars, with estimates available in the literature. On 
the other hand, for the row-shifted (2 x 2)-paramagnetic 
phase boundary at high fields, we have found a discrep- 
ancy of some 19% between our estimate and that given in 
Ref. 3, for the highest transition temperature at _ff = 6. 
Even allowing for the (rather plausible) likelihood that 
our extrapolation procedures underestimate the extent of 
the ordered phase, for PRG with the largest size avail- 
able (iV = 16) one has Tc{H = 6) = 0.651, already 11% 
below Tc{H := 6) = 0.7293(7) quoted in Ref. |3. At this 
point, such discrepancy remains unexplained. 

We have not succeeded in gathering as much infor- 
mation regarding critical properties of the row-shifted 
(2 X 2)-paramagnetic boundary line, as we did for its 



collinear-paramagnetic counterpart. However, the be- 
havior of r] illustrated in Figure [5] reminds one of the low- 
temperature behavior of the two-dimensional XY model. 
Indeed, with Tkt being the upper limit of the Kosterliz- 
Thouless critical phase, the exponent rj of the XY model 
grows smoothly and monotonically from rj = at T = 
to 1/4 at Tkt (with most of the increase confined to 
higher T: at T = 0.5Tkt, V ~ 0.05) ^. So the very 
low values of 77 found in the present case may, or may 
not, indicate the presence of incipient XF-like behavior 



along at least part of the high-field critical line. 
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